/*******************************************************************************
Figure1.do

This file identifies the 12x12 mile squares roughly containing Princeton, 
Richmond, and New York City, plots the Starbucks locations in each of these
three squares, and produces the three panels shown in Figure 1.

Throughout this file, we use the Albers equal-area ellipsoid projection to 
convert latitude-longitude coordinates into location on a plane. See file 
"Code/Setup/Continental_US_Grids.do" for more details regarding this projection.

Last updated: 11/5/2019
*******************************************************************************/

version 15
cd "C:\Plants_in_Space"
set more off
set type double
set varabbrev off

********************Locate 12x12 square containing Richmond*********************
// Use Richmond City Hall as a landmark, then use the geoinpoly command to 
// map the landmark's coordinates into the appropriate square. This square
// contains the majority of Richmond, as well as part of Chesterfield County
// to the south and Henrico County to the west.
shp2dta using "Shapefiles\Raw\Richmond\tl_2019_51_arealm.shp", ///
	database("Shapefiles\Intermediate\Richmond\Area_landmark.dta") ///
	coordinates("Shapefiles\Intermediate\Richmond\Area_landmark_coord.dta") ///
	replace
	
use "Shapefiles\Intermediate\Richmond\Area_landmark.dta", clear
keep if FULLNAME=="Richmond City Hall"
destring INTPTLAT INTPTLON, replace
drop _ID
geoinpoly INTPTLAT INTPTLON using "Shapefiles\Final\US_12_mile_square_cells_latitude_longitude.dta", unique
qui sum _ID
local id_richmond = r(mean)

********************Locate 12x12 square containing Princeton********************
// Use Princeton University as a landmark, then use the geoinpoly command to 
// map the landmark's coordinates into the appropriate square. This square 
// contains all of Princeton, as well as large parts of several surrounding
// townships.
shp2dta using "Shapefiles\Raw\Princeton\tl_2019_34_pointlm.shp", ///
	database("Shapefiles\Intermediate\Princeton\Point_landmark.dta") ///
	coordinates("Shapefiles\Intermediate\Princeton\Point_landmark_coord.dta") ///
	replace

use "Shapefiles\Intermediate\Princeton\Point_landmark.dta", clear
keep if FULLNAME=="Princeton"
assert _N==1
qui sum _ID
local id=r(mean)

use "Shapefiles\Intermediate\Princeton\Point_landmark_coord.dta", clear
keep if _ID==`id'
drop _ID
geoinpoly _Y _X using "Shapefiles\Final\US_12_mile_square_cells_latitude_longitude.dta", unique
qui sum _ID
local id_princeton = r(mean)

********************Locate 12x12 square containing Manhattan********************
// Use Manhattan's centroid as a landmark, then use the geoinpoly command to map
// the landmark's coordinates into the appropriate square. This square contains 
// most of Manhattan, a small part of Brooklyn to the southeast, and a large part 
// of New Jersey to the west.
shp2dta using "Shapefiles\Raw\NYC\tl_2019_36_pointlm.shp", ///
	database("Shapefiles\Intermediate\NYC\Point_landmark.dta") ///
	coordinates("Shapefiles\Intermediate\NYC\Point_landmark_coord.dta") ///
	replace

use "Shapefiles\Intermediate\NYC\Point_landmark.dta", clear
keep if FULLNAME=="Manhattan" // Manhattan centroid
assert _N==1
qui sum _ID
local id=r(mean)

use "Shapefiles\Intermediate\NYC\Point_landmark_coord.dta", clear
keep if _ID==`id'
drop _ID
geoinpoly _Y _X using "Shapefiles\Final\US_12_mile_square_cells_latitude_longitude.dta", unique
qui sum _ID
local id_nyc = r(mean)

*********************Keep selected squares in US shapefile**********************
// Keep the coordinates of the three squares we're interested in, project, 
// then save the cell location information for each city individually.
use "Shapefiles\Final\US_12_mile_square_cells_latitude_longitude.dta", clear
keep if _ID==`id_richmond' | _ID==`id_princeton' | _ID==`id_nyc'
geo2xy _Y _X, gen(new_Y new_X) projection(albers)
return list
local a = r(a)
local f = r(f)
local lat1 = r(lat1)
local lat2 = r(lat2)
local lat0 = r(lat0)
local lon0 = r(lon0)
drop _X _Y
rename new_X _X
rename new_Y _Y

foreach location in richmond princeton nyc {
	preserve
		keep if _ID==`id_`location''
		gen polyline=1
		save "Shapefiles\Intermediate\M12_grid_ID_`location'_projected.dta", replace
	restore
}

******************Keep Starbucks locations in selected squares******************
// Keep Starbucks locations in the squares containing Princeton, Richmond, and 
// NYC. Project the coordinates of each of these establishments, and save for 
// each city individually.
use "Data\Intermediate\NETS\NETS_HQs_at_least_5_employees_2014_cleaned.dta", replace
rename ID_12 ID
drop ID_*
keep if ID==`id_richmond' | ID==`id_princeton' | ID==`id_nyc'
keep if hq==155366107 & sic8==58120304 // HQ and SIC8 for Starbucks
save "Data\Intermediate\Starbucks\Starbucks_in_Richmond_NYC_Princeton.dta", replace

use "Data\Intermediate\Starbucks\Starbucks_in_Richmond_NYC_Princeton.dta", clear
foreach location in richmond nyc princeton {
	preserve 
		keep if ID==`id_`location''
		save "Data\Intermediate\Starbucks\Starbucks_in_`location'.dta", replace
	restore
}
geo2xy latitude longitude, gen(new_latitude new_longitude) projection(albers, `a' `f' `lat1' `lat2' `lat0' `lon0') //Units here are in meters
drop latitude longitude
rename new_latitude latitude
rename new_longitude longitude
save "Data\Intermediate\Starbucks\Starbucks_in_Richmond_NYC_Princeton_projected.dta", replace
foreach location in richmond nyc princeton {
	preserve 
		keep if ID==`id_`location''
		save "Data\Intermediate\Starbucks\Starbucks_in_`location'_projected.dta", replace
	restore
}

****************************Create map for Richmond*****************************
// For map outlines, use county lines (where Richmond City is its own county).
shp2dta using "Shapefiles\Raw\Richmond\tl_2019_us_county.shp", ///
	database("Shapefiles\Intermediate\Richmond\Counties.dta") ///
	coordinates("Shapefiles\Intermediate\Richmond\Counties_coordinates.dta") ///
	replace

use "Shapefiles\Intermediate\Richmond\Counties_coordinates.dta", clear
geo2xy _Y _X, gen(new_Y new_X) projection(albers, `a' `f' `lat1' `lat2' `lat0' `lon0')
drop _X _Y
rename new_X _X
rename new_Y _Y
save "Shapefiles\Intermediate\Richmond\Counties_coordinates_projected.dta", replace

use "Shapefiles\Intermediate\Richmond\Counties.dta", clear
destring STATEFP COUNTYFP, replace
keep if STATEFP==51
// keep Richmond city, Henrico County, and Chesterfield County
keep if COUNTYFP==760 | COUNTYFP==87 | COUNTYFP==41

grmap using "Shapefiles\Intermediate\Richmond\Counties_coordinates_projected.dta", id(_ID) osize(*0.2) ///
	point(data("Data\Intermediate\Starbucks\Starbucks_in_richmond_projected.dta") x(longitude) y(latitude) size(*0.58)) ///
	polygon(data("Shapefiles\Intermediate\M12_grid_ID_richmond_projected.dta") by(polyline)) ///
	plotregion(margin(-80 -80 -30 -20))
graph export "Plots\Figure1_Richmond.png", width(1000) replace

****************************Create map for Princeton****************************
// For map outlines, use minor civil divisions.
shp2dta using "Shapefiles\Raw\Princeton\tl_2019_34_cousub.shp", ///
	database("Shapefiles\Intermediate\Princeton\NJ_county_sub.dta") ///
	coordinates("Shapefiles\Intermediate\Princeton\NJ_county_sub_coordinates.dta") ///
	replace

use "Shapefiles\Intermediate\Princeton\NJ_county_sub_coordinates.dta", replace
geo2xy _Y _X, gen(new_Y new_X) projection(albers, `a' `f' `lat1' `lat2' `lat0' `lon0')
drop _X _Y
rename new_X _X
rename new_Y _Y
save "Shapefiles\Intermediate\Princeton\NJ_county_sub_coordinates_projected.dta", replace

use "Shapefiles\Intermediate\Princeton\NJ_county_sub.dta", clear
destring COUSUBFP, replace
// keep Hopewell township, Lawrence township, West Windsor township, Plainsboro 
// township, South Brunswick township, Franklin township, Montgomery township, 
// and Princeton
keep if COUSUBFP==33180 | COUSUBFP==39510 | COUSUBFP==80240 | COUSUBFP==59280 | ///
		COUSUBFP==68790 | COUSUBFP==24900 | COUSUBFP==47580 | COUSUBFP==60900

grmap using "Shapefiles\Intermediate\Princeton\NJ_county_sub_coordinates_projected.dta", id(_ID) osize(*0.2) ///
	point(data("Data\Intermediate\Starbucks\Starbucks_in_princeton_projected.dta") x(longitude) y(latitude) size(*0.6)) ///
	polygon(data("Shapefiles\Intermediate\M12_grid_ID_princeton_projected.dta") by(polyline)) ///
	plotregion(margin(0 0 0 -10))
graph export "Plots\Figure1_Princeton.png", width(1000) replace

*******************************Create map for NYC*******************************
// For map outlines, use major bodies of water (NY state/borough lines in many 
// cases are in the middle the water--using these bodies directly creates a more 
// clear map than using municiple boundaries). Use bodies larger than 1,000,000 
// square meters, and remove lakes, ponds, reservoirs, etc. that do not serve as 
// boundaries. 

// Water shapefile for New Jersey, Hudson County
shp2dta using "Shapefiles\Raw\NYC\tl_2019_34017_areawater.shp", ///
	database("Shapefiles\Intermediate\NYC\NJ_HudsonCounty_water.dta") ///
	coordinates("Shapefiles\Intermediate\NYC\NJ_HudsonCounty_water_coordinates.dta") ///
	replace

use "Shapefiles\Intermediate\NYC\NJ_HudsonCounty_water_coordinates.dta", clear
replace _ID=_ID+1000 // make IDs unique for when we append files for different counties together
save "Shapefiles\Intermediate\NYC\NJ_HudsonCounty_water_coordinates_edit.dta", replace
	
use "Shapefiles\Intermediate\NYC\NJ_HudsonCounty_water.dta", clear
destring AWATER, replace
keep if AWATER>1000000
drop if strpos(FULLNAME, "Lk") > 0
drop if strpos(FULLNAME, "Pond") > 0
drop if strpos(FULLNAME, "Reservoir") > 0
replace _ID=_ID+1000 // make IDs unique for when we append files for different counties together
save "Shapefiles\Intermediate\NYC\NJ_HudsonCounty_water_edit.dta", replace

// Water shapefile for New Jersey, Bergen County
shp2dta using "Shapefiles\Raw\NYC\tl_2019_34003_areawater.shp", ///
	database("Shapefiles\Intermediate\NYC\NJ_BergenCounty_water.dta") ///
	coordinates("Shapefiles\Intermediate\NYC\NJ_BergenCounty_water_coordinates.dta") ///
	replace

use "Shapefiles\Intermediate\NYC\NJ_BergenCounty_water_coordinates.dta", clear
replace _ID=_ID+2000 // make IDs unique for when we append files for different counties together
save "Shapefiles\Intermediate\NYC\NJ_BergenCounty_water_coordinates_edit.dta", replace
	
use "Shapefiles\Intermediate\NYC\NJ_BergenCounty_water.dta", clear
keep if AWATER>1000000
drop if strpos(FULLNAME, "Lk") > 0
drop if strpos(FULLNAME, "Pond") > 0
drop if strpos(FULLNAME, "Reservoir") > 0
replace _ID=_ID+2000 // make IDs unique for when we append files for different counties together
save "Shapefiles\Intermediate\NYC\NJ_BergenCounty_water_edit.dta", replace

// Water shapefile for New York, New York County (Manhattan)
shp2dta using "Shapefiles\Raw\NYC\tl_2019_36061_areawater.shp", ///
	database("Shapefiles\Intermediate\NYC\NY_water.dta") ///
	coordinates("Shapefiles\Intermediate\NYC\NY_water_coordinates.dta") ///
	replace

// Append files together so our map features all the major bodies of water from
// each of these counties, then merge individual polygons into one using mergepoly
use "Shapefiles\Intermediate\NYC\NY_water_coordinates.dta", clear
append using "Shapefiles\Intermediate\NYC\NJ_HudsonCounty_water_coordinates_edit.dta"
append using "Shapefiles\Intermediate\NYC\NJ_BergenCounty_water_coordinates_edit.dta"
geo2xy _Y _X, gen(new_Y new_X) projection(albers, `a' `f' `lat1' `lat2' `lat0' `lon0')
drop _X _Y
rename new_X _X
rename new_Y _Y
save "Shapefiles\Intermediate\NYC\NY_NJ_water_coordinates.dta", replace
	
use "Shapefiles\Intermediate\NYC\NY_water.dta", clear
destring AWATER, replace
keep if AWATER>1000000
append using "Shapefiles\Intermediate\NYC\NJ_HudsonCounty_water_edit.dta"
append using "Shapefiles\Intermediate\NYC\NJ_BergenCounty_water_edit.dta"
mergepoly using "Shapefiles\Intermediate\NYC\NY_NJ_water_coordinates.dta", ///
	coor("Shapefiles\Intermediate\NYC\NY_NJ_water_coordinates_merge.dta") replace

grmap using "Shapefiles\Intermediate\NYC\NY_NJ_water_coordinates_merge.dta", id(_ID) osize(*0.2)  ///
	point(data("Data\Intermediate\Starbucks\Starbucks_in_nyc_projected.dta") x(longitude) y(latitude)) ///
	polygon(data("Shapefiles\Intermediate\M12_grid_ID_nyc_projected.dta") by(polyline)) ///
	plotregion(margin(0 -10 -50 -50))
graph export "Plots\Figure1_NYC.png", width(1000) replace
